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We report an application of the tri-dimensional pseudo-spectral time domain algorithm, that solves 
with accuracy the nonlinear Maxwell’s equations, to predict second harmonic generation in lithium 
niobate ridge-type waveguides with high index contrast. Characteristics of the nonlinear process 
such as conversion efficiency as well as impact of the multimode character of the waveguide are 
investigated as a function of the waveguide geometry in uniformly and periodically poled medium. 


I. INTRODUCTION 


Optical nonlinear devices based on ferroelectric materials such as Lithium Niobate (LiNbOs) or Lithium Tantalate 
{LiTaOs) have a great technological potential. These cristals indeed allow the generation of any wavelength from the 
mid IR to the UV using quasi phase-matched nonlinear processes mia- This technology is usefull to realise optical 
souces or can also play a key role inside wavelength multiplexing systems for optical telecommunications. Highly 
nonlinear optical fibers are competing alternative thanks to their ability to convert wavelengths through four wave 
mixing process via the Kerr effect. However, weak third order nonlinear coefficients give weak conversion efficiency 
even over several hundred meters of propagation. To the contrary, efficient conversion efficiency are possible over 
short distance (few centimeters) in periodically poled Lithium Niobate (PPLN) due to a large second order nonlinear 
coefficient allowing fast parametric conversions. 

Particularly, ridge waveguides with high index contrast have a key role to play in the improvement of performances 
[3]. In comparison with standard waveguide fabrication techniques, such as proton exchange [4] or titanium in-diffusion 
[5] , stronger confinement can be reached with ridge waveguides along with a good overlap between fundamental modes 
of the nonlinear process thanks to a high index contrast. In addition, long term stability of the devices is improved 
especially if the nonlinear material is doped with proper elements such as magnesium to limit the photorefractive 
effect. At last, fabrication techniques based on micromachinig or etching keep the intrinsic properties of the material. 
In particular, the strong nonlinear coefficients and low optical absorption coefficient remain unchanged. However, 
high index contrasts implies very small waveguide cross sections (sub-micron square) in order to form singlemode 
waveguides. Such a strong confinment is beneficial to envision high nonlinear conversion efficiency but it implies 
tight tolerances on the fabrication process of the periodic poling for quasi-phasematching and waveguide geometrical 
uniformity. Moreover, efficient light coupling in these tiny structures is very challenging. As a consequence, an 
alternative solution is to realise waveguides with larger cross sections (from 5/im^ to 100/im^) which releases constrains 
both on the fabrication process and on light coupling. In that case, influence of the multimode character of these 
waveguides on the nonlinear process has to be studied. 

To optimize the design and performances of these waveguides, numerical modeling of waves propagation in such 
devices is essential. Many numerical methods are available for wave propagation modeling. The Finite-Difference 
Time-Domain method[8] (FDTD), the Split-Step Method[7] (SSM) and the Finite Element Method[6] (FEM) belongs 
to the time-domain methods while the Beam Propagation Method [9] (BPM), the Transfer Matrix Method (TMM) 
and the Eigen Mode Expansion method (EME) work in the frequency-domain. Among these numerical techniques, 
the EDTD (implemented in a large number of both free and commercially available sofware) is the most general 
and rigorous time-domain method. It provides solutions for a large number of guided optics configurations such 
as photonic crystal waveguides, surface plasmon waveguides, devices with high-index contrast waveguides, ring and 
disk resonators, negative index material structures, dispersive, and nonlinear materials. However, when nonlinear 
phenomena in long guiding structures are studied, the EDTD becomes prohibitively computationally intensive and 
therefore impractical. 

In order to overcome this later limitation we developed a numerical model, based on a Pseudo-Spectral Time Domain 
algorithm pT] (PSTD). The nonlinear Maxwell’s equations is solved in three spatial dimensions to fully characterise 
the electromagnetic fields of the wave injected in the ridge waveguide and the generated harmonic wave. This paper 
is organised as follows. In section [TTj the principle of the PSTD algorithm is described. Then, the numerical results of 
second harmonic generation (SHG) in uniformly poled and in periodically poled ridge waveguides are given in section 
HT| These results are compared to basic analytical calculations in order to validate the numerical model. Einally, 
modal analysis of the guided waves is performed in order to study the influence of the multimode property of the 
ridge waveguides on the SHG efficiency. 
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II. NUMERICAL MODEL 



FIG. 1: Schematic view of the modeled periodically poled lithium niobate ridge-type waveguide. The square section ridge is 
bonded to a silicon substrate with a buffer layer. 2A is the period of the poled domains and z is the propagation axis of light. 


In PSTD algorithm, Maxwell’s curl equations are calculated with discrete Fourier transforms in order to solve 
the spatial derivatives on an unstaggered, collocated grid m .This spatial differential process converges with infinite 
order of accuracy for grid-sampling densities of two or more points per wavelength m, provided that the medium 
optical properties are sampled in accordance with the Nyquist theorem. Limitations of Fast Fourier Transform (FFT), 
due to the periodic boundary conditions, are avoided by using absorbing boundary conditions formulated for Perfect 
Matched Layer (PML) in nonconductive media pTHT3] . As a consequence, this numerical method can be used to 
study various problems on larger scales, more efficiently and with a better accuracy than Finite-Difference Time- 
Domain (FDTD) methods pTl IMUTS] . In particular, because accurate modelling of nonlinear optical processes with 
the FDTD method requires extremely fine sampling to minimize numerical dispersion errors, PSTD schemes offer 
significant improvements in computational efficiency and accuracy [laiini- In the present work, 3D-PSTD algorithm 
models the SHG in periodically poled and uniformly poled Lithium Niobate (LN) ridge-type waveguides. The modeled 
device is represented schematically in Fig. 

In our numerical problem, we consider the propagation along the z axis of two transverse, electromagnetic waves 
inside the ridge-type waveguide. More specifically we consider the fundamental wave and its second harmonic (SH) 
wave in CW regime. The (x, y, z) axis correspond respectively to the X, Z, Y axis of the LN crystal {Z axis corresponds 
to the crystallographic c-axis of the LN crystal). The total volume is sampled with Ux x riy x riz points giving a 
grid-sampling density greater than two points per the shortest considered wavelength (inside the LN crystal). The 
square section waveguide is bonded to a silicon substrate with a 3 ym thick transparent buffer layer. A is the length 
of the poled domains (Fig. which gives a periodic poling of period 2A. 

Use of FFT in Maxwell’s curl equations yields time-stepping relations of the form given in reference [T7]. For 
example, the x component of the electric displacement Df and of the magnetic field Bf at the fundamental wavelength 
is expressed as: 


p \n+l _ (l—+AtFj, ( — 


( 1 ) 


where Dxyj is incremented by the partial derivative of Bzj with respect to y and Dxzj is incremented by the 
partial derivative of Byj with respect to 2 : [17]. (j, k, 1) are the indices of the cartesian coordinates at the considered 
spatial point (jAx^kAy^lAz) and n is the index of the temporal step nAt. The maximum time step to fullfill the 
stability criteria of the 3D-PSTD algorihm is At = where c is the speed of light in vacuum, is the 

boundary absorbing layer function along the v {v = x^y^z) dimension which is obtained by setting the term to 
zero inside the region of interest and to a pure imaginary part increasing linearly along the v dimension inside the 
absorbing domains m with widths of ^ sampling periods at the boundaries of the sampled volume. Fy and ^ 
correspond to the FFT and the inverse FFT along the dimension v. Other components of Df and Bf are obtained 
by circular permutation of the x, 2 : indices in Eq. Similar equations give the fields components of the SH wave. 
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The source terms that generate the fundamental wave is designed such as the transverse electromagnetic wave 
emitted by the source propagates along the z direction. They are added at the intermediate time-step n + ^, to the 
terms ^yzjYjki ^^ing the relations : 


Dyzj0= 


( 2 ) 


Sxjl^ki ^yj\]ki amplitude of the transverse components of the CW source term at the time step nAt and 

at the spatial sampling point (jAx, kAy, lAz). These components are given by: 


ifxj and ifyj define the polarization state of the electromagnetic wave emitted by the source. 5 ' 0 ,/li/C/ gives a 
gaussian beam in the (x, y) transverse plane and with the optimized three-cells normalized pattern [|, |] along the 

2 : axis in order to suppress the aliasing errors ED. Finally, the propagation of the fundamental wave in the increasing 
2 : direction is ensured by also adding the source terms to the magnetic field EH- 

In Eq. the electric field components of the fundamental and the SH waves are updated, using the full coupled 
wave equations describing the SHG process in the LN crystal, as follows: 
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where 6 / 22 , dsi and ^33 are the quadratic nonlinear coefficients of the LN crystal. These coefficients are set to zero 
outside the LN medium. The periodic poling is modelled by changing periodically the sign of the coefficient d^s (FigQ 
with the period A = Lc, where Lc is the coherence length of the SHG process. erv,w\jki is the relative permittivity 
of the sampled media along the v direction, at the spatial sampling point (/Ax, kAy, lAz) and at the wavelength 
{w = /, h). This parameter takes in account the birefringence of the LN crystal. 


III. NUMERICAL RESULTS 
A. SHG in uniformly poled LN ridge-type waveguides 

First, we use the 3D-PSTD algorithm for modeling the SHG process in an uniformly poled waveguide with a square 
section of 10 x 10 ym? . The source emits a ImW continuous (CW) gaussian beam at the fundamental wavelength 
\f which is linearly polarized along the c-axis of the LN crystal {ipf = cpxj = 0 and cpyj = 0). The gaussian beam 
profile is adjusted in order to optimize the coupling of the fundamental wave in the waveguide. The other parameters 
of the simulation are listed in table [H 

Figure depicts the intensity distribution along the ridge waveguide and at the output for both the fundamental 
and the SH waves. Dotted red lines give the position of the perfect phase matched layers (PML) at the boundaries 
of the sampled volume and dotted white lines correspond to the interfaces between the different materials (air, LN, 
buffer and silicon). Along the 30ym propagation length many features can be observed. First, over about 10ym 
length the launched beam is reshaped to form the guided beam mainly constituted of the fundamental mode. The 
oscillatory behavior of the SH beam intensity is also clearly seen with about two periods observed. At the output of 
the waveguide a wider beam size is obtained for the fundamental beam compare to the SH beam as expected from 
the wavelength difference. We would like to emphasize that the decay of the fundamental wave intensity at the end 
of the waveguide is due to the absorbing layer. Numerical simulations have also been performed for waveguides with 
5x5 and 2.5 x 2.5 ym^ square sections. 
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TABLE I: Numerical values of the simulation parameters. 


Wavelengths (in fim) 

Xf ? Xfi 

1.55, 0.775 

LN refractive indices 

'^ef ■) ’^oh-i 

2.2128, 2.1371, 2.2606, 2.1780 

LN nonlinear coefficients (in pmlV) 

^22, dsi, dss 

2.1, -4.64, -41.7 

Sampling points 

Tlx X Uy X Tlz 

128 X 128 X 256 

Spatial sampling steps 

<1 

II 

<1 

II 

<1 

Ah 

2-2ng^ 

Temporal sampling step 

At 

Ax 

8c 

Buffer layer refractive index 

Tiy 

1.501 

Silicon refractive indices 

T^Sif, TlSih 

3.50 , 3.72 

Ridge length (in fim) 

L 

30 



-20 -10 0 10 20 -10 -5 0 5 10 

z propagation axis ([im) x axis (lim) 


FIG. 2: Numerical results: (a) propagation of the fundamental wave inside a ridge-type waveguide with a square section of 
10 X 10 (b) SH wave generated in the ridge. Dotted red lines show the position of the the perfectly phase matched layers 

boundaries of the sampled volume. Dotted white lines denote the interfaces between the different materials. 


From the variation of the total power of the SH wave as a function of the propagation distance (Fig. [^, the 
coherence length of the SHG process can be estimated with a fairly good accuracy. Note that the oscillatory behavior 
is less and less perfect as the section of the waveguide is deacresed. Indeed, intensity variation of the SH wave along the 
propagation axis in the larger waveguide is almost sinusoidal in good agreement with the usual plane wave formalism 
while this is not the case for smaller waveguides even when size of the input gaussian beam is optimized. These 
discrepancies are due to the reshaping distance as well as the presence of weakly excited higher order modes. This 
latter point is confirmed by numerical simulations that show stronger intensity fluctuations of fundamental guided 
waves along the propagation axis in smaller waveguides. This suggests that injected light is coupled with higher order 
guided modes which implies that different harmonic frequency waves with different coherence lengths are generated 
simultaneously. 

Coherence lengths estimated with the PSTD simulations can be compared with the coherence lengths calculated 
with the usual definition: 


Lc 


__ 

4^{nh,eff - rif^eff) 


(6) 


where n^^eff and Tif^eff are the effective indices of the fundamental eigenmodes at the wavelengths Xh and Aj, 
respectively. For the modeled waveguides, table |H| gives the coherence lengths deduced from PSTD simulations 
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FIG. 3: For the different modeled guides, total power of the SH generated in the non poled LN ridge-type waveguides as a 
function of the propagation distance. 


and the coherence lengths calculated with Eq. where effective indices are calculated with a commercial software. 
Coherence lengths obtained from both methods are in good agreement. 

TABLE II: Comparison between the coherence lengths deduced from the PSTD simulations and calculated with Eq. 


Square section of Ridge (in ym^) 

10 X 10 

5x5 

2.5 X 2.5 

Coherence length with PSTD (ym) 

8.7 

7.7 

5.4 

LN Effective index nj^eff 

2.1340 

2.1265 

2.0975 

LN Effective index rih^ef f 

2.1773 

2.1753 

2.1677 

Coherence length with Eq. |6| {ym) 

8.9 

8.2 

6.0 


B. SHG in periodically poled LN ridge-type waveguides 

Next, the coherence lengths deduced from the above PSTD simulations (table are used to design ridge-type 
waveguides where domains are now periodically poled with a period twice the coherence length. A series of simulations 
are then performed for the three considered structures. Eigure shows the intensity distribution, in the xz and xy 
planes, of the fundamental and SH waves. Eor all square sections the transverse spatial sampling steps {Ax and 
Ay) is adjusted such as the total number of sampling points in the transverse plane is a constant. We observe that 
waveguides with smaller cross sections (Eig. and |^) show intensity fluctuations associated with shorter periods 
along propagation. Such fluctuations denote the multimode character of the guided waves eventhough the size of the 
launched gaussian beam is optimized to excite mainly the fundamental mode. As the waveguide section is reduced 
the mode beating is associated with shorter period which is in agreement with effective indices of modes that are more 
and more dissimilar. This trend is in accordance with modal properties of waveguides. As a consequence, SH waves 
also suffer from noise due to both the multimode and the phase matching conditions. Eig. shows the evolution of the 
total power of the SH wave as a function of the propagation distance. Since the total power of the launched source is 
a constant, the power density of the guided beam increases when square section of the guide decreases. Consequently, 
the nonlinear phenomena is more efficient and the total power of the harmonic wave increases. According to the basic 
theory, this gain in generated power should increases like the inverse of the section of the waveguide. This trend is 
present in the depicted results but a slightly lower gain is obtained that can partialy be attributed to the presence of 
the beam reshaping length. 

In order to estimate the efficiency of the non linear process we plot the normalized efficiencies (in %W~^cm~^) of 
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FIG. 4: Numerical results: (a,c,e) propagation of the fundamental wave through the periodically polled ridge-type waveguides 
with different square section. (b,d,f) SH wave generated in the waveguides. Dotted white lines denote the interfaces between 
the different materials. 



0 10 20 30 

z propagation axis (pm) 


FIG. 5: Total power of the SH generated waves along propagation distance in the periodically poled LN ridge-type waveguides 
of three different cross sections. 


the SHG process (Fig. as a function of the propagation length which is defined as: 


(7) 


where and Pf are the total power of the SH and fundamental waves inside the waveguides, respectively and z is 
the propagation distance. To compare with, the dotted curves correspond to the theoretical efficiencies defined for 
the first domain as : 


Vth{z) 


^^ 





( 8 ) 


where Sj are the cross section areas of the guided beam and Ak = {Tih,eff ~ '^f,eff) is the phase mismatch 
between the fundamental and the SH waves and r]max is the efficiency for perfect phase matching. Because of the 
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FIG. 6: Normalized SHG efficiencies in the periodically poled LN ridge-type waveguides. Golor dotted curves show the 
corresponding theoretical normalized efficiencies calculated with Eq. Black dotted lines symbolise for each ridge the efficiency 
level 


periodic poling of the waveguide, efficiency along the domain (with p > 1) is calculated by incrementing the 
efficiency at the end of the domain p — 1 as follows: 


llth (z) = llth {{P - l)Lc) + Vmax 


sin(^ - {p- 


1 ) 1 ) 


Akz 

2 


2 


(9) 


where an additional phase shift of — | is applied to the sinus function from one domain to the next. These curves show 
that the normalised efficiency, deduced both from the PSTD simulations and from Eq. is high at the beginning of 
the waveguide and decreases along propagation in the first domain. Then, it rapidly stabilises to an almost constant 
value. This behavior is inherent to the quasi-phase matching configuration. The reached level is well approximated 
the efficiency r]max reduced by a factor of as shown in [22]. This level is represented in Fig. [^for each ridge by the 
black dotted lines. Note that the conversion efficiency calculated from Eq. [^is slightly higher than the one deduced 
from the PSTD results since the latter method takes into account the multimode character of the waveguide. More 
importantly, these curves show that very strong conversion efficiency can be obtained when the size of the waveguide 
decreases. 


C. Modal analysis of the guided waves 

In this last section we study the influence of the multimode property of the ridge waveguides on the SHG process. We 
first exploit the numerical results given by the PSTD algorithm in periodically poled waveguides in order to evaluate 
the dominant excited eigenmode. To this purpose, a commercial software is first used to calculate the components 
of the transverse eigenmodes for each ridge waveguide at both fundamental and harmonic wavelengths. Then, for a 
given eigenmode m at the wavelength the overlap integral with the output beam obtained with the PSTD method 
is calculated using : 






PSTD ] 


PSTD 


x,w -^y,wm '^y,w 


rj* 


i) dxdy\'^ 


ff\E^^'^^\‘^dxdy P\Hwrn?dxdy 


( 10 ) 


where are the electric field components given by the PSTD algorithm and {Hx,wm,Hy,wm) are the 

magnetic field components of the eigenmode m calculated with the commercial software. As an example, Fig. |^shows 
the transverse components of the electric fields for the fundamental and the harmonic waves at the output of the 
5x5 pm? waveguide along with and the magnetic fields components of the TMqq eigenmodes. The obtained values of 
the overlap integral between the guided beam and the first mode (TMqo) at the fundamental wavelength are superior 
to 98% whatever the geometry of the waveguides. It thus shows that the guided mode is mainly composed of the first 
mode which is a favorable situation to reach high conversion efficiency. The residue of higher order modes still gives 
some intensity fluctuation due to mode beating which is noticeable for small section waveguides for which beating 
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FIG. 7: Field components of the fundamental and harmonic waves for the 5 x b iiw? ridge waveguide : (a-d) electric fields of 
the output guided waves from the PSTD algorithm and (e-h) magnetic helds of the TM fundamental eigenmodes given by the 
commercial software. Dotted white squares show the ridge contours. 


length is shorter ( see Fig. [^. Additional overlap integral calculations between the fundamental mode of the SH waves 
and the guided beam at the fundamental wavelength obtained from PSTD show that the overlap decay as the section 
of the waveguide is getting larger. For instance, the overlap for the 10 x 10 fim? ridge is 90% while it reaches 99% for 
the 2.5 X 2.5 iim^ one. Such an almost perfect matching comes from the strong similarity between the fundamental 
modes at both wavelengths as observed in Fig. This feature is an advantage of high index contrast ridge type 
waveguides. We can assert that the strong second harmonic generation efficiency reachable in narrow waveguides is 
not only due to a tight confinement but it also benefits due to a better mode matching between fundamental modes 
at both wavelengths. 

In order to complement these results, we also calculate the overlap integrals along the propagation axis. Figure 
shows the variation of the overlap integrals, between the guided waves and the TMqq eigenmodes of the waveguides, 
along the propagation axis. We note that the launched Gaussian beam shape is optimised as witnessed by an overlap 
close to 100% with the fundamental wave starting from the entrance of the waveguides. To the contrary, coupling 
with the T Mqo modes of the SH waves exhibits strong fluctuations at the early stage of the propagation and increases 
quickly up to the optimal values. This effect is especially visible in small section waveguides. We attribute the initial 
stage to the distance necessary to the reshaping and redistribution of the generated SH light among the different 
eigenmodes. 


IV. CONCLUSION 

For the first time to our knowledge, a 3D-PSTD algorithm has been implemented to solve the second harmonic 
generation process in lithium niobate ridge-type waveguides. The model has first been validated with uniformly 
poled waveguides and further used to characterise and better comprehend the physics in periodically poled ridges. 
The model is able to determine most characteristics of the SHG process such as optimum poling period, conversion 
efficiency along with accurate guided beams profile evolution. More specifically, it shows that very high conversion 
efficiency can be reached in high-index contrast ridge type waveguide even though they are not single mode. These 
performances are obtained with injection of a gaussian beam that is shaped to favour excitation of the fundamental 
mode of the guiding structure. For instance, conversion efficiency as high as 2000 %W~^ is expected in a 1cm long 
ridge waveguide of 2.5 x 2.5 jum^ square section. 
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FIG. 8: Variation of the normalised overlap integrals along the propagation axis. Overlap integrals are calculated between the 
guided light from PSTD and the first eigenmodes at the fundamental (solid curves) and harmonic (dotted curves) wavelengths 
for the 10 X 10, 5 X 5 and 2.5 x 2.5waveguides. 
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